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ABSTRACT 

The auroras on Jupiter and Saturn can be studied with a high sensitivity and resolution by 
the Hubble Space Telescope (HST) ultraviolet (UV) and far-ultraviolet (FUV) Space Tele- 
scope imaging spectrograph (STIS) and Advanced Camera for Surveys (ACS) instruments. 
We present results of automatic detection and segmentation of Jupiter's auroral emissions as 
observed by HST ACS instrument with VOronoi Image SEgmentation (VOISE). VOISE is a 
dynamic algorithm for partitioning the underlying pixel grid of an image into regions accord- 
ing to a prescribed homogeneity criterion. The algorithm consists of an iterative procedure 
that dynamically constructs a tessellation of the image plane based on a Voronoi Diagram, 
until the intensity of the underlying image within each region is classified as homogeneous. 
The computed tessellations allow the extraction of quantitative information about the auroral 
features such as mean intensity, latitudinal and longitudinal extents and length scales. These 
outputs thus represent a more automated and objective method of characterising auroral emis- 
sions than manual inspection. 
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1 INTRODUCTION 

Planetary auroral emission is a particularly useful diagnostic tool 
for the study of magnetospheric processes. In the polar regions of 
the planet, currents flow into and out of the upper atmosphere and 
energetic particles precipitate down the magnetic field lines, im- 
pacting on the molecules in the atmosphere. At Jupiter and Sat- 
urn, the dominant atmospheric species is hydrogen, which emits 
in the ultraviolet (UV) when excited by auroral electron impact. 
In order to observe this UV emission we need to use a platform 
above the Earth's absorbing atmosphere, such as the Space Tele- 
scope imaging spectrograph (STIS) instrument (until the 2004 Sat- 
urn campaign), and its successor the Advanced Camera for Surveys 
(ACS) instrument, both on board the Earth-orbiting Hubble Space 
Telescope (HST). From the study of these auroral emissions we can 
learn a great deal about magneto spheric proce sses at thes e planets 
withou t even leaving Earth orbit jClarke et alj [2002: Grod ent et al.l 
I2003bl l3). 

The auroras of the outer planets have been studied relatively 
little in comparison to the Earth's, however we have advanced our 
understanding of the environment of lupiter and Saturn. Jupiter's 
magnetosphere is the largest cavity in the solar wind flow, i.e. 
within the heliosphere, and the dynamics of it are largely domi- 
nated by the fast rotation of the planet (period ^ 10 hr). It has been 
established that the most significant component of Jupiter's auroras. 
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the main auroral oval, is driven by planetary rotation. The moon lo 
orbits deep within Jupiter's magnetosphere (orbital radius ~ 6 Rj 
where Jupiter equatorial radius Rj = 71492 km), where it is ex- 
posed to large tidal forces as it is subject to the competing gravities 
of Jupiter and the other Galilean moons: Europa, Ganymede, and 
Callisto. These tidal forces heat the interior of lo, such that it has 
become extremely volcanic, and liberates from its volcanoes of the 
order of one tonne per second of sulphur and oxygen plasma into 
the neighbourhood of its orbit. This plasma feels the fast rotation 
of Jupiter's magnetic field and is accelerated to the same rotation 
rate as the planet. This fast rotation causes it to 'diffuse' away from 
Jupiter due to the centrifugal force, and in doing so it slows down 
its rotation in an attempt to conserve angular momentum. This re- 
sults in a vast spinning disc of plasma that co-rotates with Jupiter 
in the inner region and more slowly in the outer region. Jupiter's 
magnetic field, however, attempts to keep all this plasma rotating 
at the same angular velocity. This enforcement is mediated by an 
electric current system, one component of which results in elec- 
trons bombarding Jupiter's upper atmosphere and causing the main 
oval auroras (Hill 1979; Vasvhunas 1983). 

lupiter also exhibits two other classes of auroral feature: 
the moon footprints and the polar auroras. The moon footprints 
are spots of auroral emission that are magnetically linked to the 
Galilean moons. The brightest spot is that associated with lo, and it 
also exhibits a tail, or 'wake' that traces eastward (in the direction 
of rotation) around the planet. The exact details of the processes 
that lead to the formation of the footprints are still not understood. 
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but it is thought that they are a result of the moons' interaction 
with the rotating magnetic field and the acceleration of newly cre- 
ated plasma to the spee d of the ro tating field via Alfven waves 
jGumett & Goer tz 1981: iBonfond et al. 2007). The polar auroras 
are still somewhat of an enigma, as they are the most highly vari- 
able component of Jupiter's auroras in terms of brightness. The lo- 
cation and variation of these auroral emissions make it possible that 
at least some of them may be caused by the solar wind interaction, 
althou gh this has yet to be pro ved conclusively. Interestingly, recent 
work jPallier & Prangdl200lh shows good evidence that Jupiter's 
polar cap boundary, i.e. the perimeter of the region where planetary 
field lines open into the solar wind, is coincident at times with some 
of these polar emissions. 

The preceding discussions highlight an important link be- 
tween the geometrical properties and locations of auroral features, 
and the physical processes of their origin. Previous observational 
studies, e.g. Pallier & Prange (2001), have used correlation tech- 
niques (spatial filtering) to emphasise auroral features at selected 
length scales. In this paper, we present a new automated technique 
for producing maps of the entire range of length scales and other 
properties of features in a given auroral image. 

Segmentation refers to the process of partitioning a digital im- 
age into multiple regions (sets of pixels). The goal of segmenta- 
tion is to simplify and/or change the representation of an image 
into something that is more physically meaningful and easier to 
analyse. Image segmentation is typically used to locate objects and 
boundaries (lines, curves, etc.) in images, and different approaches 
have been develop ed especially for the purpose of medical imaging 
jPham et al.ll2000h . In the astrophysical context different methods 
have also been used in image processing an d analysis. Amo ngst 
others, we mention the w atershed transform ( Platen e t alj2007h . the 
adaptive spatial binning i Cappellari & CoDin,2003;) and Bayesian 
image reconstruction dCabrera et al.ll2008h . 

In this paper we present the VOronoi Image SEgmentation 
(VOISE) algorithm, a method for segmentation of an image into 
regions based on the partitioning of the support of the image (the 
set of data values over all pixels) into a set of disjoint convex poly- 
gons. Each polygon is a Voronoi region (VR) and the boundary 
points of the collection of the VRs is a Voronoi diagram (VD). 
The VD is sometimes also known as a Dirichlet tessellation and 
the VRs are also called Dirichlet r egions, Meijering c ells, Thiessen 
polytopes, or Voronoi polygons ( IOkabeetalj|2000l) . The set of 
VRs is determined by a finite number of points (one point in- 
side each VR) called seeds, sometimes called germs or generators. 
VOISE is an iterative and self-organising algorithm (Jan tsch I980i ; 
iNicoIis & Prigoginel 19891) for automatic segmentation of an image, 
based on adaptive construction of a VD. It consists of dynamically 
subdividing and merging a network of convex polygons according 
to the information contained in the underlying image. 

In section |2l we give a detailed introduction about the VD, 
and in section[3] the VOISE algorithm is developed. The behaviour 
of auroras is frequently studied through manual inspection of two- 
dimensional image data. In section |4l we apply the VOISE algo- 
rithm to images of Jupiter's auroras from the HST ACS instrument 
with the Solar-Blind Channel (SBC). In section |5] we discuss the 
results of the segmentation. We expect the VOISE algorithm to pro- 
vide more objective and more efficient analyses than visual inspec- 
tion, especially when such analysis relies on statistics of large vol- 
umes of data, whose manual processing may be prone to human 
error. 

Additional natural byproducts of the VOISE algorithm include 
compression of data and noise-limited detection of features. These 



aspects will be explored in a future study. It is the purpose of the 
present paper to introduce the fundamentals of the technique, to- 
gether with its application in auroral physics study. 



2 VORONOI DIAGRAM 

The VD is one of a few truly interdisciplinary concepts with rel- 
evant material to be found in fields as different as anthropol- 
ogy, astronomy, e cology, physics, and urban and regional planning 
dOkabeet al][200Q) . Let us begin by considering a set of seeds 
S = {si, S2, ■ • ■ , s„}, where Si are distinct points in the plane 
R^. The VR associated with seed Si, noted TZ{si), is the set of 
points p in the plane that are nearer to Si than to any other 
seeds Sj, with respect to the Euclidean distance d, formally 

TZ{S,) ^ {p £R'\d{p,Si) <d{p,Sj), SJeS\s^}, (1) 

and its closure, noted TZ{si), is defined as 

TZ{si) = {pe^''\d{p,Si)^d{p,s,), s,eS\s,}, (2) 

The VD of a set of seeds 5*, noted V{S), is the graph or 'skeleton' 
formed by the polygonal boundaries of the VRs of all seeds in S, 
i.e. 



v{s) = fl 7^(s. 



(3) 



Since the intersection of any VRs is empty and the union of the 
closure of all VRs is the plane, the VD is called a tessellation of the 
plane. 

In the terminology of mathematical morphology, the VR 
TZisi) is called the influence zone or the dominance region of Si, 
and V{S), the skeleton by influence zones. Note that the definitions 
given m Eqs. CHUl are still valid for any space where a distance 
has been defined . For instance, in the three-dimensional Euclidean 
space, the VR are convex polyhedra and the VD consists of the 
skeleton formed by the boundaries of the VRs. 

Figure[T]illustrates the VD for a set of eight seeds in the plane. 
Two seeds are said to be neighbours if their associated VRs share a 
common edge. 

In addition to its geometrical representation, the VD can also 
be identified by an abstract representation called a Voronoi graph. 
The Voronoi graph is a planar graph in two-dimensional space. 
The graph G{S) — {V, E) represents the topological structure 
of the neighbouring VRs. The nodes (or vertices) V and lines 
(or edges) E of the graph G represent respectively the seeds and 
the unordered pairs of nodes that represent links connecting two 
seeds. In Figure [T] the neighbours of seed si consists of the set 
J\f{si) = {s2, S3, S5, S7, se}. Thus we can deduce that the pairs 
{si, S2}, {si, S3}, {si, S5}, {si, St} and {si,sb} are links of 
the graph associated to the node si. 

The data structure consisting of the neighbour list for all 
seeds 7V(5') = {Af{si),Af{s2), . . . ,Af{s„)} thus represents the 
Voronoi graph. 

The set J^{S) also characterises the Delaunay triangulation 
(DT), which is obtained by connecting all the pairs of neighbour 
seeds. The nodes and faces (triangles) of the DT correspond re- 
spectively to the polygons and nodes of the VD. The DT is called 
the dual graph representation of the primal graph of the VD. 

The DT are widely used in various applications such as finite 
element methods and interpolation schemes, owing to the fact that, 
once the set of seeds S has been defined, the DT optimises the 
geometric 'compactness' of the triangulation. In particular, the DT 
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of points constant, the process tends toward the homogeneous Pois- 
son point pr ocess and the tess ellation is called a Poisson-Voronoi 
tessellation jOkabe et all200 d). 

The set of VRs represents a 'tiling' of the image plane fl. Let 
us consider an operator f{l, Si) which associates, to each seed Si, 
a unique value depending on the region TZ{si), and the intensity 
of the image pixels within that region, i.e. the restriction 2 In^Si) 
of J to TZ{si), with values in the set {l{p),p G TZ{si)}. Then a 
tiled image 7} can be constructed by 'filling' each polygon with 
the value given by applying the operator f{l, Si), formally 



Tf{p) = f{I, Si), for Si e S and p G 7^(si 



(4) 



In order to calculate a tiled approximation of an image with a set 
of seeds S, one can for instance use statistical averages such as the 
arithmetic mean or the median for the operator /. 



Figure 1. Illustration of a planar VD (solid lines) and its dual DT (dotted 
lines) of a set 5 consisting of eight seeds (labelled 1 to 8). 



maximises the minimum angles of all its element triangles. Also 
the DT is the unique triangulation where the circumscribed circle 
of each triangle does not contain any other vertices. 

Note also that a VR may be either a bounded or an unbounded 
region. When a VR is bounded the number of edges and vertices are 
equal, when unbounded there is one vertex less than there are edges. 
In Figure [T] the regions TZ{si) and TL{s2) are bounded while all 
the others are unbounded. The polygon V consisting of the vertices 
{s4, sg, S3, S5, S7, ss}, which are the seeds of all the unbounded 
regions, is called the convex hull of the set S. 

The construction of the VD and the DT is a fundamen- 
tal problem in computational geometry, and many approaches 
have been proposed to compute them, such as the incremen- 
tal method, the divide-and-conquer method and the plane sweep 
method (Okabe et al. 2000). 

An image in this context is a real-valued function X{p) defined 
for p, pixels in the image plane: a two-dimensional lattice space 
SI C . In the rest of this paper, we restrict the set of seeds S to 
be distinct points in the image plane Q,. And accordingly, the VR 
and VD are defined as their restriction to the image plane $1. 



3 VOISE ALGORITHM 

VOISE is an iterative and self-organising algorithm, and consists of 
up to four different phases explained hereafter. 



3.1 Initialisation Ptiase 

A small number of seeds are initially randomly drawn according 
to the uniform probability distribution over the image plane Q,, and 
the corresponding VD is constructed. 

The generation of points in space according to a probability 
distribution is called a stochastic point process. For a uniform prob- 
ability distribution over a bounded region the process is the bino- 
mial point process, where for any subset of the region, the probabil- 
ity distribution for the number of points within the subset is given 
by the binomial distribution. Asymptotically, if we let the bounded 
region tend toward the infinite region while keeping the density 



3.2 Dividing ptiase 

The dividing phase is an iterative process and is the main process 
for 'evolving' the VD. First we need to define a measure of homo- 
geneity for the set of pixels of the underlying image within each 
VR of the diagram. We define the merit function x of the image X 
in the region associated to seed Si as 

max [T{p)] - min [l{p)] 
X (I, Si) = II ,^ (5) 



||x(2:,'S')|| = max 

Si£S 



Mi,s)\\ 

max [J(p)] - min [l{p)] 
peniSi) penis,) 



(6) 



The numerator of the right hand side in Eq. (|5} is a 'max-min' 
measure of the variation of the image intensity within the region 
7^(si). The numerator defined in Eq. ^ is the L-infinity norm of 
this max-min measure of the image intensity for the entire VD. The 
values of the measure of homogeneity x {T, Si) are, by definition, 
bounded in the interval [0, 1]. The minimum value of x corresponds 
to the most homogeneous VR and conversely, the maximum value 
is for the least homogeneous region. Note that the merit function x 
defined in Eqs. l[5}{6) is only one of many possibilities. Any alter- 
native merit function can easily be plugged in VOISE. 

When a cell is not homogeneous according to the chosen ho- 
mogeneity threshold Xm, i.e. 



(7) 



seeds are added within this region. The number of seeds added cor- 
responds to the number of vertices in the polygon. Each added seed 
is the barycentre of the seed of the VR with weight Ws, and a ver- 
tex of the polygon with weight Wv, corresponding to a fixed rel- 
ative distance from the seed. The most probable number of ver- 
tices/edges of a t ypical cell of a Poisson-Voronoi diagram is six 
jOkabe et alj200Ql . in other words the most probable polygon is an 
hexagon. Figure |2] illustrates the method for three iterations start- 
ing with the hexagon in the centre of the upper left panel and it- 
eratively over the innermost hexagon, adding six seeds at each it- 
eration, and for three weightings, {ws = 1/4, uii, = 3/4) in the 
upper right panel, {ws — 1/2, Wv = 1/2) in the lower let panel 
and {ws — 3/4, Wv = 1/4) in the lower right panel. In the last 
case where the seeds are added nearest the seed of the original re- 
gion, the newly created central region is very small compared to 
the newly created surrounding regions. In VOISE we set these pa- 
rameters to Ws =5/8 and — 3/8 which provides a good com- 
promise between quick convergence and comparable size of new 
regions. 
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Figure 2. Illustration of the method to add seeds for different weight- 
ings. Seeds are added to the central hexagonal region (upper left panel) 
recursively for three iterations, adding 18 seeds. The weights are Ws = 
1/4, = 3/ A (upper right panel), Ws = 1/2, Wy = 1/2 (lower left 
panel) and Ws = 3/4, = 1/4 (lower left panel). 

This iterative and deterministic algorithm of adding new seeds 
ensures the regular growth of regions such that self-similar struc- 
ture across different length scales is emphasised (fractals with re- 
cursive definition such as the Koch flake are classic examples of 
self-similarity). 

The threshold Xm is a dynamic and self-organising parameter 
estimated for each new iteration of a dividing phase using quantile 
statistics on the sample Xi of the measure of homogeneity, i.e. 

X. = {x(X,sO,s, GS}. (8) 

Xm is defined as the value of the homogeneity measure for which 
a fraction po/lOO of the polygons are homogeneous, i.e. such that 
x(lj ^i) < Xm- Thus defining the probability Prob, that a given 
polygon has its value x^X-^^i) below the threshold Xm, we can 
write the following 

Prob(x(J,sO<Xm) = ^ (9) 

Simultaneously, the probability that a given polygon has its value 
x(l) ^0 exceeding the threshold Xm is given by 

Prob(x(X,s,:)^Xm) = l-^ (10) 

The value Xm is obtained at each iteration as the inverse of the 
cumulative histogram of the sample Xi for the value pD (see middle 
panel in Figure[5l(. Since the cumulative histogram is a monotonic 
increasing function, smaller values for pd give smaller value for 
the threshold Xm and conversely. 

Note that as the subdividing process takes place, the norm 
||x(lj'S')|| is becoming smaller, and the measure of homogeneity 
is becoming more stringent. The set of homogeneity measures will 
accumulate towards the superior value sup(x) = 1- The dividing 
phase is iterated until all VRs are classified as homogeneous, or 
the size of the polygon does not allow the further addition of new 
seeds. This is characterised by enforcing a minimum distance dm 
between two seeds. 

For a large value po of the percentile (70-80 %), the thresh- 
old Xm increases, and thus the dividing scheme will only apply to 
the (100 — Pd) % most inhomogeneous regions according to the 
merit function x- With the merit function defined as in Eq. (O, a 



large value of pD will only allow the dividing process to lock onto 
the most contrasted regions and will not necessarily be able to de- 
tect low-contrasted or noisy features of the image. Conversely, for 
a small value pd of the percentile (60-70 %), the dividing process 
will be more systematic and uniform over the entire image but will 
not account for the most contrasted features of the image. However, 
it is important to emphasise that the VOISE algorithm is not lim- 
ited to the present merit function and any other merit function that 
would help characterise a feature, e.g. texture, noise properties, can 
be plugged in. 

The minimum distance dm between seeds is the parameter 
which limits the length scale of the features one wishes to detect. 
This parameter should be chosen according to the size of the small- 
est resolvable features. 

3.3 Merging phase 

At this stage, any VR is either homogeneous with respect to the 
prescribed threshold Xm, or its shape does not allow to add any new 
seed as described in the previous section and fulfil the minimum 
distance requirement dm between neighbouring seeds. 

The merging phase is also an iterative process. It consists of 
removing unnecessary 'fine' regions where adjacent VR have very 
similar characteristics, e.g. the averaged intensity fj,i of the image 
for the region associated to seed Si as calculated by Eq. {4} with 
a specified operator /, in the following the median of the set of 
values. 

A seed Si is defined as redundant and ought to be removed 
when the following conditions are fulfilled: 

(i) The region TZ{si) is homogeneous with respect to a chosen 
homogeneity threshold Xm, i-e. 

x(X,s.)<Xm (11) 

The value for the threshold is estimated at each iteration from the 
cumulative histogram of the merit function sample Xi from a pre- 
scribed PM-th percentile, in the same way as for po in section IT2l 

Prob(x(J,sO<Xm) = ^ (12) 

(ii) amongst its neighbouring seeds Sj £ N{si) that are embed- 
ded in homogeneous region, the relative difference between the av- 
eraged values fj,i and fj,j is less than a prescribed threshold A/i, 
formally 

l^i - Mil < ^mIm^I- (13) 

Two such regions are said to be identical with respect to the aver- 
aged intensity ^. 

(iii) the total length C of the edges shared with non homogeneous 
neighbours normalised to the perimeter V of the region does not 
exceed a prescribed threshold A7i, formally 

C/V < AH. (14) 

At each iteration, all seeds are checked for redundancy and 
removed when fulfilling the three conditions above. After a set of 
seeds has been removed, some regions have grown to absorb the 
seeds that have been removed. In such regions, the merit function 
X can only increase (and may become non homogeneous), and the 
averaged intensity is changing. Thus during the merging process, 
the three conditions become more difficult to fulfil for any seeds 
and the process stops by itself ('self-regulating'). 

The percentile of the merging phase, pm, does not need to 
be identical to the percentile of the dividing phase, pD, but should 
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be larger in order to relax the homogeneity criteria and conversely 
smaller in order to strengthen the homogeneity criteria (see Fig- 
ure [Sjl. 

The merging phase does not disturb the overall representation 
of the image but provides a better approximation with fewer seeds. 
When a seed is removed, the neighbouring polygons reorganise 
themselves to fill the region associated to the removed seed. We 
noted that the polygons tend to organise according to the geometri- 
cal shapes in the image. The final organisation is made independent 
of the initialisation phase by using the merging phase. 

Our tests on the image presented in the next section show that 
different initial distributions of seeds do not significantly change 
the final VD. Typically the numbers of polygons, i.e. the number 
of seeds in the final VD is affected by less than 4-5 %. The typical 
difference in area of the VRs for different realisations is less than 
10-30 pixels^. And the typical difference in intensity is less that 
5-10 %. Here we computed the difference in area of the VRs be- 
tween two realisations as follows. In a given realisation Ri of the 
VOISE algorithm, each polygon TZis}) with seed s} has an area 
A} . To compute the difference in polygon area between Ri and a 
second realisation R2, we determine for each seed s} from Ri, the 
polygon TZis]) from R2, with seed s'j and area which contains 
this seed, i.e. such that s\ G TZ{s^). The difference in area is then 
Ai . 
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3.4 Regularisation phase 

The tessellation obtained from the divide-and-merge phase can 
be 'regularised' through an iterative relaxation process. This is 
achieved by computing iteratively a new VD where each seed Si 
is replaced by the mass centroid (centre-of-mass) of the polygon 
TZ{si) with a density function p 



E ppip) 
E p(p) 

PS-R.(S,) 



(15) 



When the position of each seed coincides with the centre-of- 
mass of the polygonal region, the partitioning is called a centroidal 
Voronoi diagram (CVD) tessellation and the method is known as 
Lloyd's method. The CVD corresponds to an optimally regular and 
uniform tessellation with a minimum-energy configuration, in the 
sense that it minimises the norm of the second-order inertial mo- 
ment of each polygon. Asym ptotically, the C VD approaches a uni- 
form hexagonal-like lattice iIDu etal.lll999l) for constant density 
function. 

The density funct ion p has an im portant role in the regular- 
isation as it is seen in IDu et alj ( Il999h . Let us consider the case 
where the image itself is used as density function, i.e. p{p) is the 
value of the image data in the pixel at position p. After the divid- 
ing and merging phases, the polygons are distributed such that they 
minimise the merit function x defined in Eqs. l|5}{6]l, and the associ- 
ated density function should be, to a good approximation, constant 
within each polygon. 

Note that the regularisation process is controlled in two man- 
ners. The iteration process stops by itself when the relative distance 
between the centre-of-mass and the seed of any region, i.e. 



-s,| for all s,eS, 



(16) 
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is less than one pixel. One pixel is enough since the positions of the 
centres-of-mass are rounded off to an integer number in order to be 



Figure 3. Sub-image extracted from the polar projection of Jupiter's UV 
auroral emission computed from an HST ACS/SBC image taken on Febru- 
ary 21, 2007 at 16:03:58 UT. The size of the image is 256 X 256 pixels 
and the axis units are expressed in Jupiter's equatorial radius (71492 km). 
The origin corresponds to the projection of the northern rotational pole of 
Jupiter and the colour-coded intensity is in kR. 



in the image plane fl. In addition, and if required, it is possible to 
specify a maximum number of iterations of the relaxation process. 

The regularisation scheme could also be used as a precondi- 
tioning to the first iteration(s) of the dividing phase. Following the 
initialisation phase where a limited number of seeds are drawn ran- 
domly, a regularisation phase with the image as density function 
tends to move the seeds to positions where the image data tends to 
form local 'peaks'. 



3.5 Implementation 

Note that due to the dynamic nature of the algorithm where seeds 
and added and removed, it would be quite inefficient to construct 
the entire VD every time the seed configuration is modified, and a 
more efficient approach is to use an incremental algorithm which 
operates only on a given seed and its neighbours (see Figure |2]l- 
In addition the VRs are needed for the points in the image plane, 
which consists of a discrete and well organised set of points. We 
have chosen to use an incremental algorithm for discrete VD us- 
ing a structure similar to the one described in lSegueira & Preteuxl 
l ll997l) . which is designed for seeds located on the image plane f2, 
i.e. at discrete positions. 

VOISE requires a set of parameters to be initialised. These 
consist of pu and dm for the dividing phase, pM, A/i and AH for 
the merging phase, and an optional maximum number of iterations 
for the regularisation phase, corresponding to partial regularisation. 
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4 APPLICATION TO HST IMAGE 

Jupiter's UV auroras have been observed for many years but re- 
cently high-resolution images using the HST STIS and ACS in- 
struments have revealed highly dynamic auroral emissions th at are 
controlled to a significan t degree by the solar wind ( Clarke et_alj 
l2002l;lNichols et alj2008l) . 

Figure [5] shows an image of Jupiter's UV auroral emis- 
sion in the Northern hemisphere taken on February 21, 2007 at 
16:03:58 UT and projected onto a northern polar view. 

The raw data were collected with the SBC channel of the ACS 
through the F125LP long-pass filter, with short wavelength cutoff 
of 125 nm. This filter mostly excludes the H Lyman-a band. The 
instrument consists of a 1024 x 1024 Multi Anode MicroChannel 
Array detector with an average scale of ~ 0.032 arc sec pixel"^, 
such that the overall field of view is 35 x 31 arc sec^. During each 
100s exposure the blurring introduced by planetary rotation of any 
co-rotating features is ~ 1 deg at the Central Meridian Longitude 
(CML). The raw images are reduced: corrected for geometric dis- 
tortion and scaled to 0.0250 arc sec pixeP^, flat-fielded and dark- 
count calibrated using the latest calibration files available from 
the Space Telescope Science Institute. The images are then con- 
verted from counts pixeP^ to kR of H2 and Lyman-a emission 
(where 1 kR represents a photon source flux of 10^ cm~^ s~^ ra- 
diating into 47r steradians) usi ng the conversion f actor 1 kR — 
1.473 • 10"^ counts s"^ pixel ^Nichols et all2008l) . 

The image in Figure [3] is a reduced polar projection com- 
puted using the meta data provided with the original image and 
the NA SA Navigatio n and Ancillary Information Facility SPICE 
system jActorjl 19961) . and assuming the light emission is from an 
infinitesimally thin shell located 200 km above the 1 Bar level. The 
1 Bar level corresponds to a oblate spheroid of eccentricity 0.3 and 
semi-major axis 71492 km. 

The CML of Jupiter at the time of the observation is 
166.77 deg. The orientation of the image is such that the s-axis 
corresponds to the CML meridian, and the dimensions along the 
axes are in units of Jupiter's radius (R,j = 71492 km). 

A non-continuous auroral oval is clearly identified, as well as 
several patches of emission within associated to plasma injection 
and particle precipitation. The footprint of lo is seen in the upper 
right part of the image at coordinates (O.GRj, 0.15Rj). The foot- 
print of the magnetospheric cusp is seen at (0.4Rj, — 0.05R,j). 

Figure |4] illustrates the result for the three phases: dividing, 
merging and regularisation of the VOISE algorithm applied to the 
HST image presented in Figure [3] The 'tiled' images are con- 
structed by filling each polygon with the median value of the pixels 
that lie within each polygon, and the VD consists of the black solid 
lines. The VD was initialised with twelve random seeds. 

The upper image shows the result of the dividing phase. The 
minimum distance between two seeds was set to = 7 pixels^ 
and the homogeneity of the VD by a requirement of a percentile 
Pd ~ 85 %. The percentile po = 85 % is directly read from the 
ordinate of the cumulative histogram of the sample Xi in the mid- 
dle panel of Figure |5] and the corresponding Xm is read from the 
abscissa. The legend gives the iteration number and the number 
of seeds for the corresponding histogram. The dividing phase con- 
verged with the prescribed parameters in 27 iterations and about 
950 seeds, as seen in the top panel of Figure [5] Initially the cu- 
mulative histogram has the signature of a uniform probability dis- 
tribution, and as the dividing process takes place it converges to a 
probability distribution where most of the regions (po — 85 %) 
have their homogeneity criteria x < Xm ''^ 0.25. 
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Figure 4. Illustration of the VOISE segmentation algorithm on the image 
shown in Figure [3] Top panel: result of the dividing phase; middle panel: 
result of the merging phase; bottom panel: result of the regularisation phase. 
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The middle image shows the result of the merging phase. The 
threshold for relative similarity is = 20 %, the ratio of non- 
homogeneous boundaries is set to ATI = 30 % and the homogene- 
ity of the VD is set to a prescribed percentile pM ~ 50 %. As for 
the dividing phase, the percentile po is directly read in the ordinate 
of the cumulative histogram of the merit function x in the middle 
panel of Figure |5] and the corresponding Xm is read in the abscis- 
sas. The merging phase consists of seven iterations after which the 
process stops by itself, as no seeds can fulfil the redundancy condi- 
tions described. 

The lower image shows the result of the regularisation phase. 
The tessellation results from two relaxation iterations, i.e. replacing 
the seed of each polygon by the centre-of-mass of the associated 
region twice. Our experience shows that one or two iterations are 
in general enough, due to the rounding off of the coordinates of the 
centre-of-mass, given by Eq. J15b . to integer values corresponding 
to points of the image plane Q. 

It is worth noting that the VOISE algorithm based on a sta- 
tistical distribution of seeds evolves dynamically according to self- 
organising of the information in the image as represented by the 
changing homogeneity threshold Xm ■ 



5 DISCUSSION 

Note that any quantity or measure could be calculated for each VR, 
based on either the properties of the polygons, the statistical prop- 
erties of the underlying image or a combination of both. 

Many analytical properties of the VR have been determined 
assuming a Poisson-Voronoi tessellation. These include the prob- 
ability distribution of the perimeter segments, of the angles at the 
vertices, the polygon area, the polyg on to be n-sided, to name a few 
dOkabe et ani200(]| : lHilhorstll2008r) . Any deviations of the statis- 



tics of a property of the VD from the probability distribution of 
the Poisson-Voronoi tessellation is an indication of seed distribu- 
tion related to some rea l point pattern, or clustering as defined in 
iKaufman & Rousseeuwl ( l2005h . 

Such 'tiled' images can thus be used to extract quantitative 
information about structures present in the image, but also poten- 
tially to clean an image from unwanted parts such as contamination 
or defects. 

Figure |6] shows the length scale calculated as the square root 
of the number of pixels in each Voronoi polygon times the edge of 
a pixel in units of 100 km at the planet. This map of length scales 
allows post-processing such as selection of points or regions of the 
image according to a prescribed rule. 

Figure|7]shows the spatial distribution of the seeds of the VD 
of the middle panel of Figure |4] in Cartesian and polar coordinate 
systems. Let us consider a group of neighbouring seeds, i.e. where 
each polygon in the group shares at least one edge with another 
polygon. Such group of seeds is defined as a cluster if the polygons 
in the group are similar with respect to a given property, or combi- 
nation of properties (Kaufman & Rousseeuw 2005). Such proper- 
ties could be, for example, purely geometrical (e.g. area, perimeter, 
length scale), or related to data values from the underlying image 
inside polygons (e.g. median intensity, mean intensity). It is im- 
portant to point out that geometrical properties need not be used, 
although they would be a natural choice for detection of auroral fea- 
tures. Four such clusters with similar median intensity and length 
scales are identifiable and annotated in the polar coordinate system: 
(i) the footprint of lo, (ii) the dusk side segment of the main auroral 
oval, (iii) the dawn side segment of this oval, and (iv) a structure 
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Figure 5. Top panel: the number of seed.s, card(S'), as function of iteration; 
middle panel: the cumulative histogram of the sample Xi associated to the 
homogeneity merit function x, as defined in Eqs. ^5](6). for a selection of 
iterations of dividing and merging phases. These iterations are labelled in 
the legend together with the number of seeds; bottom panel: the threshold 
Xm derived from the cumulative histogram with the prescribed po = 85 % 
and pm = 50 % as function of iteration. Solid lines with cross markers 
refer to dividing phase data while solid lines with diamond markers refer to 
merging phase data. 



associated to cusp emission. The number of polygons per unit area 
could be considered a crude representation of the 'density of infor- 
mation' for the image. 

Our experiments with mapping polar emissions have fre- 
quently detected arc-like auroral features inside the main oval, but 
running almost parallel to it. These will be the subject of another 
study. However they are simil ar to the auroral signatur es of the po- 
lar cap boundary discussed bv lPallier & Prangj ( l200ll) . 

It is interesting to note that the clustering 'compactness' of 
the seeds is dependent on the coordinate systems. For instance the 
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Figure 6. Length scale estimated from the geometry of the Voronoi poly- 
gons obtained after regularisation (bottom panel in Figure |4). The length 
scale within a VR is defined as the square root of the area of the polygon. 
The length scale is colour-coded in logarithmic scale. 



footprint of lo in the upper right comer appears more clustered in 
polar coordinates. The scattered and relatively low number of seeds 
allows the high compression of information (with loss) achieved 
by the VOISE tessellation. The original size of the image shown 
in Figure [3] is 256 x 256 pixels with pixel intensity coded over 
4 bytes, thus requiring 256^ • 4 bytes. The tessellation contains 
877 points. Coding the median intensity over 4 bytes and each co- 
ordinate over 2 bytes requires 877 • (4+2+2) bytes, leading to a 
compression factor of 37; alternatively coding the coordinates over 
4 bytes reduces the compression factor to 25. Quantitative esti- 
mates of parameters or features such as discontinuities in the au- 
roral oval, bifurcation of auroral arcs, polar cap position, footprint 
of moons, could easily be obtained by fitting the scattered relevant 
seeds to an analytical model. Some of these features are indicated 
in Figure [T] The automatic detection of the limb of the planet to 
locate accurately and objectively the centre of the planet disc is an- 
other potential for application of the VOISE algorithm. 

A common concern that runs through a wide range of dis- 
ciplines is the examination of the spatial occurrence of a partic- 
ular phenomenon, or point pa tterns, such as clustering analysis 
jKaufman & Rouss eeuw|[2005h . The analysis of point patterns of 
the seeds (or spatial occurrence) can be done to divide the VD 
into component parts. For instance it was suggested that individual 
points may be considered to belong to one of five types: isolated 
points, members of a curvilinear structure, members of a cluster 
with an empty interior, or mem bers of ei ther the boundary or the 
interior of a non-empty cluster l lOkabe e t al. 2000). The edges of 
the polygons here aid this classification. 

Another application is multi-temporal or multi-source image 
analysis and especially the detection of event defined as a structure 
modification between two images. This can be achieved by com- 
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Figure 7. Spatial distribution of the seeds for the segmentation obtained 
the merging phase {middle panel in Figure|4) expressed in Cartesian {upper 
panel) and in polar {lower panel) coordinate systems. 



paring the VD of two images generated by VOISE according to the 
following definition. 

Let us consider two sets of seeds 5*1 and S2 and their associ- 
ated VD. Two seeds Si £ Si and Sj £ S2 are equivalent seeds 
when the seed Si belongs to TZisj), and simultaneously the seed 
Sj belongs to TZ{si). The regions associated to seeds Si and Sj are 
identical or equivalent when the averaged intensities fii and fij of 
the two images for the region associated to seed Si and Sj are simi- 
lar with respect to an imposed threshold, while an event is detected 
between the two images if the averaged intensities are not similar. 
The application to analysis of auroral features and their temporal 
variability is evident, and will be pursued in future studies. 



6 CONCLUSION 

We have presented VOISE, a self-organising dynamic algorithm 
for the automatic segmentation of an image based on adaptive con- 
struction of a VD according to some information contained in the 
underlying image. 

The VOISE algorithm has been applied to the detection of UV 
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auroral emission regions on Jupiter as observed by tlie ACS/SBC 
instrument on board the HST. 

In a following paper, we plan to do a statistical survey of 
the length scales, morphologies and variability of the auroras for 
an extensive set of HST ACS/SBC images of Jupiter, and look 
for scale-invariant characteristics, a property of dynamic systems 
known as self-organised criticality. At present, sequences of typi- 
cally about 18 images (covering time intervals about or less than 
2hrs) are available for 55 different days during the year 2007, 
providing many different views of both the northern and southern 
UV auroras. The geometrical nature of auroral emissions is known 
to be linked with their physical origin. Bright concentrated 'arcs' 
are the result of strong electric currents flowing between the iono- 
sphere and magnetosphere, while more extended, diffuse emissions 
arise from scattering of the magnetospheric particles onto trajec- 
tories which intercept the planetary atmosphere. Since the macro- 
scopic behaviour of the solar wind displays spatial and temporal 
scal e-invariance which are characteristic of self-organ ised critical- 
ity jChapman etalll 19981 : Ichapman & Watkinsll200lh . one would 
expect similar macroscopic behaviour of the auroral activity to ex- 
hibit related spatial and temporal scale-invariance. 

The VOISE algorithm could easily be extended to three- 
dimensional data (discrete volume composed of voxels) to perform 
automatic volume segmentation. 

The applicability of the VOISE algorithm is not limited to seg- 
mentation of auroral images but can in practice be used for any 
types of image data with a suitable homogeneity merit function. 
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